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Ocean  acoustic  tomography  is  a  method  for  measuring  the  sound  speed  structure  of 
a  volume  of  ocean  indirectly  using  the  transmission  of  sound  between  instruments 
deployed  within  the  area  of  interest.  Since  acoustic  tomography  is  based  on  the 
measurement  of  travel  time  of  transmitted  signals,  the  data  are  very  sensitive  to  the 
relative  displacement  of  the  instruments.  This  paper  describes  a  navigation  system 
which  continuously  tracks  ocean  acoustic  tomography  instruments.  By  this  system, 
the  contaminating  effect  of  mooring  motion  on  the  measurement  of  travel  times  of 
acoustic  signals  is  removed. 
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Introduction 


Ocean  acoustic  tomography  is  a  method  for  measuring  the  sound  speed  structure  of  a  volume 
of  ocean  indirectly  using  the  transmission  of  sound  between  instruments  deployed  within  the  area 
of  interest.  The  method  can  provide  useful  space  and  time  resolution,  and,  it  requires  less  ship 
time  than  that  presently  required  to  deploy  and  recover  conventional  oceanographic  instruments. 
Acoustic  tomography  takes  advantage  of  the  facts  that  low  frequency  sound  can  travel  long  distances 
in  the  ocean,  and,  that  travel  time  is  a  function  of  temperature  and  current  velocity.  Travel  times 
of  acoustic  pulses  can  be  interpreted,  using  inverse  theory,  to  provide  information  about  the  oceanic 
circulation  complete  with  its  random  motions. 


Since  acoustic  tomography  is  based  on  the  measurement  of  travel  times  of  transmitted  signals, 
the  data  are  very  sensitive  to  the  relative  displacement  of  the  instruments.  Given  a  typical  sound 
speed  of  1500m/sec,  15  m  of  error  in  the  estimation  of  the  length  of  an  acoustic  ray  adds  10  msec 
of  travel  time  error.  This  is  important  when  compared  with  40  msec,  the  expected  order  of  travel 
time  changes  due  to  the  intermediate,  or  oceanic  mesoscale  held.  In  other  words,  mooring  motion 
can  introduce  a  variability  about  the  travel  times  which  can  swamp  the  mesoscale  variations.  The 
structure  of  a  typical  acoustic  mooring  is  shown  in  Fig.l. 


Fig.l  Structure  of  an  acoustic  mooring  at  rest  depth.  Note  that  a  horisontal  shift,  as  a  result  of 
mooring  motion,  is  accompanied  by  a  deepening  of  the  instrument. 
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The  mooring  motion  is  monitored  by  means  of  an  acoustic  navigation  system.  This  system 
consists  of  a  micro-processor  controlled  transceiver  which  is  mounted  as  part  of  the  sensor  package, 
and  three  recoverable  transponders  anchored  on  the  bottom,  approximately  2000  m  from  the 
mooring  anchor.  At  pre  determined  intervals  the  transceiver  simultaneously  interrogates  the  three 
transponders  and  measures  the  round-trip  travel  times.  These  data  are  recorded  along  with  the 
time  and  date  so  that  after  processing,  the  path  of  the  mooring  may  be  reconstructed.  A  scheme 
of  the  geometry  is  shown  in  Fig. 2.  A  similar  system  has  been  described  previously  by  Nowak  and 
Mealy. (1) 


vu*  -  s 

no*r 


Fig. 2  Mooring  navigation  scheme 


Inferences  in  the  oceanographic  field  by  acoustic  tomography  techniques  ate  questionable  unless 
the  recorded  arrival  times,  which  are  used  in  tomography  to  infer  oceanic  structure,  are  known  with 
respect  to  the  position  of  the  instrument.  The  arrival  times  between  arrays  of  sources  and  receivers 
are  separate  and  distinct  from  the  arrival  times  between  the  three  transponders  and  the  transceiver 
which  are  used  to  track  the  motion  of  a  mooring.  It  is  the  latter  type  of  arrival  times  which  are  the 
subject  of  this  report. 

This  technical  note  reports  on  the  work  which  was  done  as  a  contribution  to  a  French/US 
ocean  acoustic  tomography  experiment  which  was  conducted  in  the  Gulf  Stream  by  scientists  from 
IFREMER/WHOI  in  1987  (the  results  of  which  are  yet  to  be  published).  This  note  details  an 
analytical  expression  which  may  be  used,  for  speed  and  simplicity,  as  a  close  approximation  to  the 
exact  or  ray  trace  solution  for  the  purpose  of  monitoring  mooring  motion. 


Method 


Once  the  path  of  a  ray  has  been  traced  from  the  transceiver  to  a  transponder  it  is  possible 
to  calculate  the  travel  time,  by  integrating  along  the  ray  path.  A  sound  speed  profile,  which  is 
representative  of  the  Gulf  Stream,  is  shown  in  Fig. 3.  This  profile  was  input  to  a  ray  trace  program 
to  find  arrival  times  for  horizontal  ranges  between  500  m  and  4500  m. 
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Fig, 3  Typical  Gulf  Stream  sound  speed  profile 

Arrival  time  (r)  may  be  expressed  in  terms  of  its  dependent  variables,  the  depth  of  the 
transceiver  (zs)  and  the  horizontal  range  between  the  transceiver  and  the  transponder  (r),  thus 
r(r,  z5).  The  ray  trace  program  was  run  several  times  for  selected  source  depths,  ranging  from 
1000  m  to  2000  m.  Here  depth  is  assumed  to  increase  downwards  with  sero  depth  occurring  at 
the  surface.  The  arrival  times  obtained  are  said  to  represent  the  exact  solutions.  It  is  possible 
to  construct  a  curve  which  relates  arrival  lime  to  horizontal  range  for  a  particular  source  depth. 
Examples  of  three  such  curves  are  drawn  Fig. 4. 


Fig  4  Arrival  time  is  plotted  against  horizontal  range  for  various  source  depths. 
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Obviously,  an  infinite  number  of  range  solutions  is  possible  for  any  arrival  time  if  depth  infor¬ 
mation  is  unavailable.  However,  if  three  arrival  times  are  known  for  a  single  given  source  depth,  it 
is  possible  to  interpolate  to  determine  horisontai  range.  A  faster  method  of  determining  horizontal 
range,  given  arrival  time  and  source  depth,  is  to  derive  an  analytical  expression  which  approximates 
the  exact  solution. 

The  methodology  used  for  determining  horisontai  range  was  to  make  a  comparison  between 
approximations  obtained  from  three  different  analytical  expressions  with  results  obtained  from  the 
ray  trace  solution.  These  expressions  for  arrival  time,  approximate  the  sound  speed  profile,  and 
or,  the  exact  ray  path.  It  is  convenient  in  this  section  of  the  report  to  refer  to  Table  1,  in  which 
each  expression  is  placed  according  to  these  approximations.  Expressions  1  and  2  were  found  to 
b«  unsatisfactory  approximations  to  the  exact  solution  for  the  purposes  of  mooring  motion.  Their 
derivation  is  given  in  Annexes  1  and  2  respectively,  for  completeness  and  for  use  in  other  possible 
applications.  All  parameters  shown  in  expressions  1  and  2  are  simil&rily  defined  in  Annexes  1  and 
2. 

Expression  3  proved  to  be  a  close  approximation  to  the  exact  solution  and  its  derivation  and 
validation  follows. 

Table  1  Analytical  expressions  which  were  tested  as  approximations  to  the  ray  trace  solution 


Sound  Speed  Profile  Approximation  Exact  Solution 


Ray  Path 


Approximation 


1  r 


- ^4 _ 

C.-fcfc 


3  r  =  T.v/l  +  (£)J 


Exact  Solution  2  r  =  ^  In  Ray  Trace 


Approximation  to  Ray  Trace  Solution 

Consider  a  ray  travelling  from  S  to  R  through  the  ocean  volume  in  which  the  speed  of  sound 
C,  varies  with  depth,  i.  The  ray  leaves  the  source  at  a  depth,  25  and  an  angle,  <*,,  where  a,  is  the 
initial  angle  as  measured  with  respect  to  the  vertical.  The  relevant  geometry  is  given  in  Fig. 5 


Fig.5  Sketch  of  how  a  ray  trajectory  may  be  approximated  by  a  straight  ray  path 
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The  integral  travel  time,  between  the  source  and  the  receiver  in  the  vertical  direction,  (T, )  is 
found  by  integrating  the  sound  speed  profile  from  the  source  depth  to  the  receiver  depth.  Assuming 
a  straight  ray  path,  the  ray  reaches  R,  the  depth  of  the  transponder,  after  an  integral  travel  time 
which  is  given  by: 

rji  _  fn  da 

~  Jz,  C(x) 

The  travel  time  is  described  by 

T  -  f  jU- 
T  ~  JzmC(x) 

where  da  is  the  distance  measured  along  the  ray,  and  the  integral  is  taken  along  a  straight  ray  path. 
Since 

*  =  ssfej 

the  arrival  time  may  be  written  with  respect  to  the  distance  along  the  ray  SR 

r  -  4  r*  ~d4~ 

rox(m)  J a  s  Cf*) 

Now  substituting  for  T„ 
r  = 

=  v^W 

where  r  =  h  tan(a«) 
and  h  —  zs  -  zr 
thus  r  =  T.  yj  1  +  ( j*)2 

Simpson’s  composite  algorithm  is  called  to  approximate  the  integral  travel  time,  T.  using  a 
number  of  sub-intervals,  m.  Annex  3  gives  details  of  this  numerical  integration  technique  and  values 
of  7i  for  source  depths  between  1000  m  and  3000  m.  The  water  depth  is  assumed  to  be  5000  m-  Fig 
6  presents  the  results  graphically. 


Fig6  Graph  showing  the  near  linear  relation  between  the  integral  travel  time,  T.  and  source  depth, 

*5- 
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The  analytical  expression  r  =  T„  yT+Tj)2  closely  approximates  the  exact  solution  for  sc  uice 
depths  between  1000  m  and  3000  m  and  horizontal  ranges  between  1000  m  and  5000  m.  A  comparison 
was  made  between  arrival  times  obtained  with  this  expression  and  those  obtained  by  using  the  ray 
trace  program.  The  largest  difference  in  arrival  times  between  the  two  solutions  is  found  at  5000 
m  range  (0.6  msec).  This  represents  a  difference  in  range  of  approximately  0.9  m.  In  all  cases  the 
differences  in  travel  time  are  small.  It  is  predictable  that  as  the  ranges  are  extended  further,  a 
straight  ray  trajectory  will  become  an  increasingly  worse  approximation  to  the  exact  solution. 

The  mooring  tracking  system  is  able  to  tolerate  the  size  of  the  differences  observed,  which  are 
in  the  order  of  tenths  of  milliseconds.  In  fact,  the  system  is  more  likely  to  be  limited  by  propagation 
loss  out  to  5000  m,  than  by  the  use  of  an  approximation  to  the  ray  trace  solution. 


Arrival  Time  Expressed  in  Terms  of  Instrument  Position 

Once  the  arrival  time  has  been  closely  approximated  by  an  analytical  expression,  it  only  remains 
necessary  to  introduce  the  other  information  which  is  available  to  the  navigation  system,  that  is, 
the  known  positions  of  the  transponders  to  arrive  at  a  solution  for  the  mooring  position.  Although 
it  is  impossible  to  assign  positions  “a  priori”  to  transponders,  once  deployed, the  positions  of  the 
transponders(r*,  yty  zx)  are  determined  by  survey.  For  our  purposes  we  will  assume  that  the  anchor 
is  at  the  origin.  Additional  assumptions  are  that  the  ocean  bottom  is  flat  and  5000  m  below  the 
surface. 

The  arrival  times  are  input  into  an  iterative  program  which  gives  the  position  of  the  transceiver 
in  terms  of  its  cartesian  coordinates.  The  program  solves  for  x,y  and  z  by  using  the  following 
equation  where  i  =  1,  2,  3. 


The  iterative  routine  finds  a  sero  of  a  system  of  three  functions  in  x,  y ,  z  by  a  modification  of  the 
Powell  hybrid  method(2).  This  routine  is  further  explained  in  Annex  4.  which  provides  a  numerical 
solution  for  three  non-linear  equations  relating  arrival  time  to  instrument  position.  Tv  is  calculated 
for  any  depth  between  1200  m  and  2000  m  by  means  of  a  cubic  spline  interpolation.  A  listing  of 
the  output  from  the  cubic  spline  interpolation  is  given  in  Annex  5.  The  FORTRAN  source  code  for 
the  mooring  tracking  program  is  given  in  Annex  6  and  a  summary  of  the  inputs  and  the  processing 
steps  which  are  required  to  achieve  a  solution  is  given  in  Fig7.  Numerical  Algorithms  Group(NAG) 
maths  library  routines(3)  are  called  to  perform  both  the  iteration  and  the  interpolation. 

Curve  Fitting  for  the  Integral  Travel  Time  Term 

The  problem  of  fitting  a  curve  for  the  integral  travel  time  term  is  solved  by  approximating 
the  values  for  the  integral  travel  time  between  transceiver  and  transponder  to  a  function,  /(x).  In 
addition,  it  would  be  advantageous  if  the  navigation  system  was  able  to  make  a  close  approximation 
of  the  instrument  position  if  depth  information  was  not  known.  A  close  approximation  was  made 
by  fitting  the  travel  time  data  to  a  least  squares  cubic  spline  curve.  A  least  squares  cubic  spline  is 
used  in  preference  to  a  least  squares  polynomial  fit  since  there  is  only  one  independent  variable. 

Under  normal  routine  operation  of  the  mooring  tracking  system,  depth  information  is  available 
from  the  pressure  sensor  which  is  mounted  as  part  of  the  instrument  package  and  the  navigation 
system  will  call  the  cubic  spline  routine  to  evaluate  Tv.  However,  if  the  pressure  sensor  fails  so  that 
depth  information  is  lost,  a  back-up  system  is  available  to  replace  the  spline  interpolation. 
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In  this  case,  since  the  required  function  is  almost  linear  and  extremely  well-behaved,  one  can 
adopt  an  expression  of  the  form  A  -f  Bz,  in  order  to  determine  T,  -  A  single  linear  equation  is  valid 
for  a  distance  of  approximately  200  m.  For  example,  assume  that  a  mooring  which  has  a  rest  depth 
of  1500  m  does  not  move  deeper  than  1700  m,  then  Tm  may  be  approximated  to  a  linear  fit  for 
1700  >  zs  >  1500  a iT,  =  3.31010  —  6.69  x  10-4  zs.  Further  examples  of  linear  fits  between 
certain  possible  intervals  are  given  in  Table  2. 


Fig7  Scheme  of  the  inputs  and  processing  required  to  arrive  at  a  solution  for  transceiver  position. 
Solid  lines  represent  main  path  for  which  source  depth  is  accurately  known.  Stippled  lines  indicate 
back-up  path  for  which  source  depth  is  approximately  known. 
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Rest  Depth  (m) 


A 


B 


1200 

3.3136 

-6.705 

X 

10- 

1300 

3.3402 

-6.895 

X 

10- 

1400 

3.3109 

-6.695 

X 

10- 

1500 

3.3101 

-6.690 

X 

io- 

1600 

3.3101 

-6.890 

X 

10* 

1700 

3.3092 

-6.685 

X 

io- 

1800 

3.3074 

-6.675 

X 

10* 

Table  2.  Coefficients  obtained  by  fitting  the  integral  travel  time  term,  Tv  to  the  source  depth,  z$, 
linearly 

The  advantage  of  fittin  T,  linearly,  in  the  absence  of  depth  information,  is  that  the  iteration 
routine  is  able  to  solve  for  u  a  third  unknown,  in  each  of  the  three  equations.  A  close  estimate 
of  the  instrument  position  is  possible  because  the  coefficients  A  and  B  change  slowly  with  depth 
because  the  sound  speed  profile  is  close  to  linear  at  these  depths. 

To  reiterate,  with  normal  navigation  system  operation,  that  is,  complete  with  instrument  depth 
information,  curve  fitting  for  T,  against  zs  is  best  done  with  a  least  squares  cubic  spline  routine. 
If  depth  information  is  known  only  to  within  200  m  accuracy  it  is  necessary  to  use  a  routine  which 
relies  on  a  close  linear  fit  for  T,  versus  2s- 


Results 

A  simulated  survey  problem  with  transponders  at  Bi  (1000, 800, n)  B 3  (1100,-900,0)  and 
By  (-1200,  -200,0)  was  designed  to  test  the  mooring  tracking  system.  Fig.  8  shows  the  instru¬ 
ment  in  its  unperturbed  position  at  Af  (0,0,1200).  The  ray  trace  solution  was  used  to  derive 
arrival  times  between  each  of  the  3  transponders  and  the  transceiver.  Accurate  positioning  of 
the  instrument  was  obtained  when  the  instrument  was  without  any  mooring  motion,  as  shown  by 
x(-0.055),y(0.055), 2(1199.34) 


Fig  8.  Transponder  positions  with  respect  to  the  mooring  anchor. 

The  moorings  were  then  perturbed  with  respect  to  lean  angle  and  lean  direction.  For  both 
maximum  generality  and  simplicity  the  moorings  were  assumed  to  be  rigid  and  the  geometry  follows 
the  lead  given  by  Cornuelle  (4).  Fig.  9  shows  ;ae  instrument  position  changes  as  a  result  of  mooring 
lean. 


FLOAT 


Fig.  9.  Mooring  lean  geometry  Az  =  /  sin  cos  Ay  =  /  sin0sint/>,  Az 
from  Cornuelle(4) 


Fig.  10  presents  the  results  of  this  simulation.  Ranges  are  from  400  m  to  2000  m.  The  ray 
tracing  for  this  work  was  done  by  specifying  the  starting  point  and  initial  direction  of  the  ray  and 
treating  it  as  an  initial  value  problem.  With  this  particular  ray  tracing  program  there  was  no 
control  over  the  point  of  emergence  of  any  particular  ray.  To  overcome  this  problem,  the  position 
of  the  receivers  was  adjusted  to  ensure  that  an  arrival  time  was  determined  for  all  lean  angles  and 
directions.  Arrival  times  and  ranges  were  evaluated  to  milliseconds  and  centimetres  respectively. 
The  system  accuracy  suffered  if  the  measurements  were  not  made  to  this  precision.  The  extremely 
small  percentage  errors  show  that  the  tracking  system  was  able  to  accurately  keep  track  of  the 
mooring  motion  even  when  the  mooring  was  subjected  to  significant  lean. 


Fig  10a  Results  of  simulation,  percentage  error  versus  lean  direction  for  x  coordinate. 
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Fig.  10b  Results  of  simulation,  percentage  error  versus  lean  direction  for  y  coordinate. 


Fig.  10c  Results  of  simulation,  percentage  error  versus  lean  direction  for  z  coordinate. 


Acknowledgements 

This  work  has  been  realised  during  my  stay  as  a  Defence  Fellow  in  the  Dept,  of  Ocean  Studies, 
Centre  Oceanologique  de  Bretagne,  IFREMER,  during  which  time  Fabienne  Gaillard  and  Yves 
esaubies  gave  much  assistance  and  kindness.  The  support  provided  by  officers  from  Director 
General  Training  and  Education  Policy,  the  administrators  of  the  Defence  Fellowship  Scheme,  is 
also  gratefully  acknowledged. 


References 

1.  Nowak,  R.T.  and  Mealy,  S.F.,  “Mooring  Motion  Monitoring  Acoustic  Navigation  System"  ,  Proc. 
IEEE  Conf.  Oceans’81,  793  -  796. 

2.  Powell,  M.J.D.,  A  hybrid  method  for  nonlinear  algebraic  equations.  In  ‘Numerical  Methods  for 
Nonlinear  Algebraic  Equations’, Ed.  Rabinowits,  P.,  Gordon  and  Breach, 1970. 

3  NAG  Fortran  Library  Manual  (NAGFLIB:  1982/1784:Mkll:Jan  84)  Numerical  Algorithms  Group 
Ltd  Manual  Vol.  1. 

4.  Cornuelle,  B.,  “Inverse  methods  and  results  from  the  1981  ocean  acoustic  tomography  experi¬ 
ment".  Ph  D  thesis.MIT/WHOI  359pp. 

5.  Clay,  C.S.  and  Medwin,  H.,  “Acoustical  Oceanography:  Principles  and  Applications".  1977. 
Wiley  Interscience. 


12 


\ 


Annex  1:  An  analytical  expression  which  approximates  the  ray  trace  solution  by  assuming 
approximations  to  both  the  sound  speed  profile  and  the  ray  trajectory. 


Expression  1  of  Table  1  takes  an  approximation  for  the  sound  speed  profile  with  a  straight  line 
approximation  for  the  ray  path.  To  obtain  an  expression  for  arrival  time,  an  incremental  step  in  the 
ray  trajectory,  ds,  is  taken,  ds  =  co^‘a  j 

A  linear  approximation  for  the  sound  speed  gradient  may  be  written  in  the  form 

C,  =  C„  +  bz 


where  C,  is  the  sound  velocity  as  a  function  of  depth,  6  is  the  sound  speed  gradient,  and  C„  is  a 
constant. 

Arrival  time,  r  =  ^;j  /* 

where  h.  is  the  distance  between  the  transceiver  and  the  transponder. 


Thus  r  : 


4ru»(aJ  C„ 


l"(g.) 


bh  =  C„(  1  _ 

horizontal  range,  r,  is  thus 
r  =  ^  tan(o,)  1  -  e 


The  required  inputs  for  this  expression  are  an  initial  launch  angle  and  arrival  time,  the  expression 
calculates  the  horizontal  shift  of  the  transceiver  for  a  given  depth.  Fig  A1  shows  the  differences  in 
range  estimates  between  an  expression  which  approximates  both  the  sound  speed  profile  and  the 
ray  path,  and  the  exact  solution.  Two  source  depths  are  drawn.  Range  differences  are  large  for 
all  initial  launch  angles.  The  expression  is  therefore  unsuitable  to  use  for  the  purposes  of  mooring 
motion. 
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Annex  2:  An  analytical  erpression  which  Appiojumates  the  ray  trace  solution  by  assuming  an 
approximation  to  the  sound  speed  profile. 

Expression  2  of  Table  1  is  a  ray  trace  approximation  which  uses  an  exact  ray  path  and  a  linear 
approximation  for  the  sound  speed  profile.  The  sound  speed  approximation  takes  the  form 

Ct  =  C„  +  bz 

where,  C„  is  a  constant  and  6  is  the  sound  speed  gradient 

*=  4c£n- 

Ct  is  the  speed  of  sound  at  depth  z.  The  expression  for  arrival  time  is  taken  from  Clay  and  Med  win 
(5).  Fig  A21  shows  the  circular  ray  path  for  a  linear  dependence  of  sound  speed  on  depth. 


The  expression  requires  a  range  and  arbitrary  a,  as  input,  to  determine  arrival  time.  In  order 
to  find  eigenrays  for  the  derived  arrival  time,  the  angles  a,  and  Qj  are  put  into  a  second  expression 
which  gives  the  horisontal  range,  r,  assuming  the  true  ray  trajectory. 

r  =  [co»(e»,)  - 

where 


and 


C,  is  the  speed  of  sound  at  the  transceiver  depth. 


c 


Fig  A2.1  Circular  ray  path  for  a  linear  dependence  of  sound  speed  on  depth. 

The  expression  may  then  be  used  to  construct  arrival  times  versus  range  curves  for  various 
transceiver  depths.  However,  the  range  scales  are  such  that  it  is  difficult  to  obtain  sufficient  resolution 
to  determine  the  difference  in  range  between  the  approximation  and  the  trace  from  these  curves. 
The  difference  in  range  is  best  shown  if  a  direct  comparison  is  made  between  ranges  and  arrival 
times  obtained  with  the  expression  with  those  obtained  after  interpolation  of  the  ray  trace  solution. 
Table  A2  presents  examples  of  this  comparison. 

Source  Depth  =  1500  m 
Speed  of  Sound  —  1493.1  m/s 


Analytical 

Ray 

Range 

Expression 

Trace 

Difference  (m) 

2047.828  m 

2037.5  m 

10.33 

2  669  s 

2.669  s 

3629.403  m 

3618.75  m 

10.66 

3.318  s 

3.318  s 

4172.547  m 

4162.5  m 

10.05 

3.584  s 

3.584  s 

4726.363  m 

4718.75  m 

7.61 

3.870  s 

3.870  s 
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Source  Depth  =  1000  m 
Speed  of  Sound  =  1492.8  m/s 


Analytical 

Ray 

Range 

Expression 

Trace 

Difference  (m) 

712.954  m 

625.00  m 

87.95 

2.674  9 

2.674  9 

1530.816  m 

1487.50  m 

43.32 

2.819  a 

2.819  9 

3085.750  m 

3056.25  m 

29.50 

3.325  s 

3.325  9 

4691.562  m 

4662.50  m 

29.06 

4.057  9 

4.057  9 

As  expected,  the  difference  in  range  is  greater  for  shallower  transceiver  depths.  Clearly,  the 
differences  in  horizontal  range  between  the  approximation  and  the  exact  solution  are  too  large  for 
this  method  to  be  considered  as  a  possible  option  for  use  in  tracking  mooring  motion. 
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Annex  3:  A  numerical  integration  technique  used  to  evaluate  the  integral  travel  time  term,  T». 


Several  numerical  methods  are  available  for  integrating  the  sound  speed  profile  between  z$  and 
zn\  the  method  which  is  used  here  is  Simpson’s  composite  algorithm.  The  procedure  is  to  divide  the 
interval  into  n  sub-intervals  and  use  Simpson’s  rule  on  each  pair  of  consecutive  sub-intervals.  Since 
each  application  of  Simpson’s  rule  requires  2  intervals,  n  must  be  an  even  integer,  that  is,  n  =  2m 
for  some  integer  m. 

With 


and 

z$  =  2 0  <  •••  <  22 m  =  ZR 

where 

Xj  —  Xi)  +  jh.  for  each  j  =  0,  1, 2m 

f  f(x)dx  =  J2  f  f(x)dx 

J*n  j-i  j 

^  /(*)<<*  =  £  +  4/1*2,.,)  +  /(*2,)J  -  6 

where  6  is  an  error  term. 

Table  A3  gives  values  for  the  integral  travel  time  between  transceiver  and  transponder  assuming 
a  straight  ray  path  and  a  water  depth  of  5000  m. 

Table  A3 


Source  Depth  (zs)  No.  af  sub-intervals  (m)  Integral  travel  time  (T.) 


1000 

200 

2.6417 

1200 

190 

2.5077 

1400 

180 

2.3736 

1600 

170 

2.2397 

1800 

160 

2.1059 

2000 

150 

1.9724 

2200 

140 

1.8390 

2400 

130 

1.7059 

2600 

120 

1.5731 

2800 

110 

1.4405 

3000 

100 

1.3081 

Annex  4:  A  numerical  solution  for  three  non-linear  equations  relating  arrival  time  to  instrument 
position. 

From  the  analytical  expression  which  closely  approximates  the  ray  trace  solution 


w 


1  + 


and  distance  formula 

r,2  =  (s  -  i,)2  +  (y  -  Vi)2 
substitute  for  arrival  time 


for  t  =  1,2,3 

So,  for  instrument  position  X ,Y,Z 


X=(x-  x2f  -  (z„  -  2S)2(£)2  +  (*„  -  zs)2  +  (V  -  V,)2 

r  =  (y-  y3)2  -  (m  -  *s)2(£)2  +  (*r  -  ‘s)2  +  (*  -  *,)2 


z  = 


■'’V1  + 


n 


An  iterative  routine  finds  a  sero  of  a  system  of  3  equations  in  3  variables  by  a  modification 
of  the  Powell  hybrid  method.  It  requires  as  input,  initial  estimates  for  X  and  Y  which  are  within 
approximately  500  m  of  the  true  instrument  position. 

Since  the  term  Tv  is  depth  dependent,  a  close  estimate  for  zs  is  essential.  This  routine  chooses 
the  correction  at  each  step  as  a  convex  combination  of  the  Newton  and  scaled  gradient  directions. 
Under  reasonable  conditions,  this  guarantees  global  convergence  for  starting  points  far  from  the 
solution  and  a  fast  rate  of  convergence.  The  Jacobian  is  updated  by  the  rank  -1  method  of  Broyden. 
At  the  starting  point  the  Jacobian  is  approximated  by  forward  differences,  but  these  are  not  used 
again  unless  the  rank-1  method  fails  to  produce  satisfactory  progress 


Annex  5:  Cubic  spline  fit  to  arbitrary  data. 

Input  Data 

Number  of  data  points  =  31  Number  of  intervals  =  4 
Unit  weighting  factors 


R 

ABSCISSA  X(R) 

ORDINATE  Y(R) 

1 

1200 

2.5077 

2 

1220 

2.4943 

3 

1240 

2.4809 

4 

1280 

2.4674 

5 

1280 

2.4540 

a 

1300 

2.4406 

7 

1320 

2.4272 

8 

1340 

2.4138 

9 

1360 

2.4004 

10 

1380 

2.3870 

11 

1400 

2.3736 

12 

1420 

2.3602 

12 

1440 

2.3468 

13 

1480 

2.3334 

14 

1480 

2.3200 

15 

1500 

2.3066 

18 

1520 

2.2932 

17 

1540 

2.2798 

18 

1560 

2.2664 

19 

1580 

2.2530 

20 

1600 

2.2397 

21 

1820 

2.2263 

22 

1640 

2.2129 

23 

1660 

2.1995 

24 

1680 

2.1861 

25 

1700 

2.1728 

28 

1720 

2.1594 

27 

1740 

2.1480 

28 

1760 

2.1328 

29 

1780 

2.1193 

30 

1780 

2.1193 

31 

1800 

2.1059 

Results 


J 

KNOT  K(J  +  2) 

B  SPLINE  COEFF  C(J) 

1 

2.5077 

2 

1200 

2.4853 

3 

1300 

2.4183 

4 

1300 

2.3066 

5 

1700 

2.1950 

6 

1800 

2.1282 

7 

2.1059 
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Cubic  spline  approximation  and  residuals 


ABSCISSA 

1200 

1220 

1240 

1260 

1280 

1300 

1320 

1340 

1360 

1380 

1400 

1420 

1440 

1460 

1480 

1500 

1520 

1540 

1580 

1580 

1580 

1600 

1620 

1640 

1680 

1680 

1700 

1720 

1740 

1760 

1780 

1800 


APPROXIMATION 

RESIDUAL 

2.5077 

■  22489E  — 6 

2.4943 

45562E-6 

2.4809 

-55100E— 7 

2.4674 

-  88447E-7 

2.4540 

26823E-6 

2.4406 

27963E  — 6 

2.4272 

90167E-7 

2.4138 

-  19607E-6 

2.4004 

-  48508E-6 

2.3870 

31715E-6 

2.3738 

—  42864E— 6 

2.3602 

—  42884E-6 

2.3468 

21135E-6 

2.3334 

31860E— 6 

2.3200 

-.  12862E-7 

2.3066 

-  68905E-6 

2.2932 

—  .62783E— 6 

2.2798 

20535E-6 

2.2664 

— .  16684E  — 6 

2.2664 

-.16684E-8 

2.2530 

27828E— 6 

2.2397 

56338E-6 

2.2263 

71112E— 6 

2.2129 

74419E— 6 

2.1995 

-.31475E-6 

2.1861 

—  44302E-6 

2.1728 

—  61795E— 6 

2.1594 

73370E-6 

2.1460 

-  35173E-6 

2.1326 

49676E-6 

2.1193 

19922E-5 

2.1059 

-.10024E-6 
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Annex  6:  A  listing  of  the  mooring  tracking  program.  The  source  code  is  written  in  FORTRAN 
for  an  IBM  PC. 

C  Program  TRACK. FOR 

C 

C 

C  This  program  finds  a  sero  of  a  system  of  3  non-linear 

C  functions  in  3  variables  by  a  modification  of  the  Powell 

C  hybrid  method. 

C 

C  Subroutine  SPLINE  is  called  to  find  Tv  for  a  given  source 
C  depth  betwen  1200m  and  2000m. 

C 

C  N  -  integer  specifying  the  number  of  equations. 

C  X  -  real  array  of  DIM(N),  it  contains  the  point  at  which 

C  the  functions  are  to  be  evaluated.  Before  entry,  X(j) 

C  must  be  set  to  a  guess  at  the  jth  component  of  the 

C  solution  ( j  =  1,2, ...n).  On  exit  j  contains  the 

C  final  estimate  of  the  solution  vector. 

C  FVEC  -  real  array  of  DIM(N),  it  must  contain  the 

C  value  of  the  ith  fn.  evaluated  at  the  point  X. 

C  IFLAG  -  integer,  if  program  termination  is  required 

C  it  is  set  negative. 

C  TOL  -  real,  XTOL  must  specify  the  accuracy  in  X  to 

C  which  the  solution  is  required,  the 

C  recommended  value  is  the  square  root  of  the 

C  machine  precision. 

C  WA  -  real  array  of  DIM(LWA),  used  as  workspace. 

C  LWA  -  integer,  LWA  =  1  N  *  (3  *  N  +  13) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

INTEGER  J.N.NOUT 

COMMON  /ITER/  Tl,T2,T3,Tv,Xl,X2,X3,Yl,Y2,Y3,Rl 
DOUBLE  PRECISION  WA(33) 

CHARACTER*!  LF.CR 

DOUBLE  PRECISION  DSQRT,Tl,T2,T3,Tv,Xl,X2,X3,Yl,Y2,Y3 
DOUBLE  PRECISION  R1,FNORM,TOL,F05ABF,X02AAF,X(3),FVEC(3) 
EXTERNAL  FCN 
DATA  NOUT  /16 / 

OPEN  (16, FILE  =  ‘TRACK  RES’,  STATUS  =  ‘NEW’) 

WRITE(NOUT,999) 

LF  =  CHAR(10) 

CR  =  CHAR(  13) 

C 

WRITE(  NOUT, 999) 


WRITE(V)TNPUT  THE  3  ARRIVAL  TIMES  (sec)’,  LF 
READ  (*,*)  T1,T2,T3 

WRITE(*,*)TNPUT  THE  OCEAN  DEPTH  (m)’,LF 
READ(V)  R1 

WRITE)', ’('INPUT  THE  X  AND  Y  COORDINATES  OF  THE-’ 

WRITE)', *)‘3  TRANSPONDERS,  ASSUME  THE  ANCHOR;’ 

WRITE)', ')‘1S  AT  THE  ORIGIN  e  x.  1000,-1400;’ 

WRITE(\*)‘  -900,1100;’ 

WRITE(*,*)‘  -1400, -1600’, LF 
READ)',')  Xl,Yl,X2,Y2,X3,Y3 
C 

C  THE  FOLLOWING  STARTING  VALUES  GIVE  A  ROUGH  SOLUTION 

C 

WRITE(  V)TNPUT  THE  LAST  KNOWN  X  COORDINATE;’ 

WRITE)  V)‘OF  THE  TRANSCEIVER  (m)’,  LF 
WRITE(*,*)TNPUT  THE  LAST  KNOWN  Y  COORDINATE;’ 

WRITE(  V)‘OF  THE  TRANSCEIVER  (m)’,  LF 

WRITE)', '(‘INPUT  THE  DEPTH  OF  THE  TRANSCEIVER  (m)’,  LF 

READ)*,')  X(3) 

C 

CALL  SPLINE(X(3),Y) 

Tv  =  Y 
N  =  3 

TOL  =  DSQRT(X02AAF(0.00) 

1FAIL  =  0 
LWA  =  33 

CALL  C05NBF(FCN,N,X,FVEC,TOL,WA,LWA,IFAIL) 

FNORM  =  F05ABF(FVEC,3) 

WRITE(NOUT,998)  FNORM,IFAIL,(X(J).  J  =  1,3) 

WRITE  (NOUT.996)  Tv,Rl,Tl,T2,T3,Xl,Yl,X2,Y2,X3,Y3 
STOP 

999  FORMAT(4(lX/),36H  TRACK-MOOR  EXAMPLE  PROGRAM  RESULTS  /IX) 

998  FORMAT(5X,31H  FINAL  L2  NORM  OF  THE  RESIDUALS, E12.4//5X, 

1  15H  EXIT  PARAMETER, I10//5X.27H  FINAL  APPROXIMATE  SOLUTION// 

1  (5X,  3E12.8) 

996  FORMAT(/6H  Tv  =,E20.5/6H  R1  =  ,E20.0/6H  TI  =,E20.5/6H 

1  T2  =,E20  5/6H  T3  =,E20.5/9H  XI, Y1  =,2E20  5/9H 
1  X2,Y2  =,2E20.5/9H  X3,Y3  =,2E20.5/) 

END 

C 

SUBROUTINE  SPLINE(XARG.FIT) 

C 

C  This  routine  computes  a  weighted  least  squares 
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C  approximation  to  an  arbitrary  set  of  data  points 

C  by  a  cubic  spline  with  knots  prescribed  by  the  user 

C 

IMPLICIT  DOUBLE  PRECISION  (A-K.O-Z) 

INTEOER  NOUT,M,NCAP,NCAP2,NCAP3,NCAP7,J 
INTEGER  IFAIL.R.NIN 

DOUBLE  PRECISION  SS,XARG,FIT,X(41),Y(41),W(41) 
DOUBLE  PRECISION  C(41),WORK(41),WORK2(4,41),K(4I) 
DATA  NCUT  /16/.NIN  /25/ 

OPEN(25,FILE  =  ‘SPLINE2.INP’, STATUS  =  ‘OLD’) 

IFAIL  =  0 

C  M  is  the  number  of  data  points 

C  NCAP  is  the  number  of  pre-defined  intervals 

C  The  knots  are  in  the  array  K 

C 

READ(NIN, 99997)  NCAP 
NCAP2  =  NCAP  +  2 
NCAP3  =  NCAP  +  3 
NCAP7  =  NCAP  +  7 
IF  (NCAP  EQ.  1)  GO  TO  40 
READ(NIN, 99995)  X(R),Y(R) 

W(R)  =  1.0 
80  CONTINUE 

IFAIL  =  1 

CALL  E02BBF(NCAP7,K,C,XARG, FIT, IFAIL) 

IF  (IFAIL  .NE.  0)  GO  TO  300 
GO  TO  320 

300  WRITE  (NOUT, 99979)  R.XARG 

320  CONTINUE 

RETURN 
C 

99997  FORMAT(I4) 

99996  FORMAT(E21. 1 1) 

99995  FORMAT(2E21.11) 

99979  FORMAT(lH,E20.5,23H  ARGUMENT  OUTSIDE  RANGE) 

END 

C 

SUBROUTINE  FCN(N,X,FVEC,IFLAG) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

INTEGER  IFLAG.N 

DOUBLE  PRECISION  FVEC(N),X(N) 

COMMON  /ITER/  Tl,T2,T3,Tv,Xl,X2,X3,YI,Y2,Y3,Rl 
C 
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X(l)  =  X 
X(2)  =  Y 
X(3)  =  Z 

FVEC(l)  =  (Tv*DSQRT(l+((X(l)-Xl)**2+(X(2)-Yl) 
1  **2)(Rl-X(3))**2))-Tl 

FVEC(2)=(X(1)-X2)**2-(R1-X(3))**2*((T2/Tv)**2) 

1  +  (Rl-X(3))**2+(X(2)-Y2)**2 
FVEC(3)=(X(2)-Y3)**2-(Rl-X(3))**2*((T3/Tv)**2) 

1  +  (Rl-X(3))**2+(X(1)-X(3)**2 

RETURN 

END 
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